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Abstract 

Recently the regular condit ional distributions of max-infinitely divisible processes were derived by 
IDombrv and E vi-Minkol 1201 ill and although these conditional distributions have complicated closed 
forms. iDombrv et all [201l| introduce an algorithm to get conditional realizations of Brown-Resnick 
processes. I n this paper we derive the regular condition al distributions of the max-stable process in- 
troduced bv lSchlatheJ [20021 ] and adapt the framework of lDombrv et al.l [20 111 ] to this specific process. 
We test the methods on simulated data and give an application to extreme temperatures in Switzer- 
land. Results show that the proposed sampling scheme provide accurate conditional simulations and 
can handle real-sized problems. 

Some key words: Conditional simulation, Max-stable process, Regular conditional distribution, 
Schlather process, Spatial extreme, Temperature. 
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Max-stable processes play a major role in the statistical modelling of spatial extremes since they arise 
as an extension of the extreme value theory to the infinite dimensional case, i.e., extremes of stochastic 
processes. Motivated by the need of suitable models for the areal modelling extremes of environmental 
processes such as heat wav es or rainfall, t he las t dec ade has seen many ad vances to develop a geostatistics 
of extremes. For instance Coolev et al. [2006j and Naveau et al.l 120091 develop variogram like tools to 
assess the spatial dependence of extremes while Padoan et al.l 2010j propose an inferential procedure 
based on the maximum pairwise likelihood estimator to fit max-stable processes to extreme spatial data 
sets. A broad over view of the different approaches available for the statistical modelling of spatial 
extremes is given bv lDavison et al. I [201l| . 

Although conditional s imulations of Gaussian based stochastic processes are known for a long time 
Chiles and Delfinerl . [l999] . conditional si mulation of max-stable processe s is a long standing problem and 



a first attempt dates back to the works of Davis and Resnick 1993 . 1989j . Recently this problem enjoyed 
renewed interest with the pioneer work of Wang and Stoev [2011 1 who give an efficient procedure to get 
condi tional simulations of max-stable processes with discrete spectral measures. Dombry and Evi-Minkol 
[2011 derive the (regular) conditional distribution of max-infinitely di visible processes and an algorithm 
to get conditional realizations of Brown-Resnick processes is given by Dombry et al.l [201 1 1 . 



Although theoretical expressions exist for the general framework of max-infinitely divisible processes, 
there are only few models where these expressions are tractable and where efficient simulation algorithms 
are possible. The aim of this paper i s to get conditio nal simulations for the extremal Gaussian process 



also known as the Schlather process [Schlatherl . 12002 



= V2n max Q max{0, £j(x)}, 



X G 
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where {Ci}i>i are the points of a Poisson process on (0, 00) with intensity measure dA(C) = £ 2 d£ and £j 
are independent copies of a standard Gaussian process with correlation function p. The processes £j are 
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assumed to be independent of the points {Ci}j>i- It is well known that the process O is a max-stabl e 
process with unit Frechet margins whose finite dimensional distribution functions are [Schlatherl 12002] 

max{0, e(xj)} 



Pr[Z(x 1 ) < zi, . . . , Z(xk) < Zk] 
for all k G N, xx,...,Xk G R d , zi, ■ ■ 



cxp 



-E 



'2ir max 
j=i,...,k 



,Zk > 0, and that its extremal coefficient function is 
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An important property of the Schlather process is that it cannot reach independence at long distances 
since 6(h) — > 1 + \ZV2 as II ^11 00 ■ Recall that indepe ndence is reached when 9(h) = 2 and that when 
9(h) — s- 2 as \\h\\ oo the max-stable process is mixing Kabluchko and Schlather , 20 ld| . 

Given a study region X C M. d , our goal is to sample from Z | {Z(x±) = z\, ■ ■ ■ , Z(xk) = Zk} for some 
(xi,...,x k ) G X k , k G N, and (zi,...,z fc ) G (0,oo) fc . InSectiond] we derive the (regular) conditiona l 



distribution of the Schlather process and recall the sampling scheme introduced bv lDombrv et al.l [2011 1 . 
Section [3] analyzes the performance of this algorithm on simulated data. The paper ends with an 
application on extreme temperatures in Switzerland followed by a brief discussion. 



2 Conditional simulation of Schlather processes 



This section aims at deriving the regular conditional distribution of Schlather processes. For our purposes, 
instead of working with the spectral characterization (TTJ), it is more convenient to use the following 
representation 

Z(x) = v27rmaxCi£i(a;), (2) 



i>l 



where {Q}i>i and are defined as in (TT|). These two spectral characterizations are equivalent since 
Z corresponds to pointwise maxima over an infinite number of random functions and consequently the 
negative parts of the Gaussian processes Si have no contribution. More formally for x G X, we have 



E 



{v / 2¥Ce(x)>0} 



Pt{Z(x) < 0] = exp 
and for k G N, (xi, . . . , Xk) G X k and zi, . . . , Zk > 
Pr[Z(xi) < Zi,.. .,Z(xk) < z k ] = cxp <^ - 
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exp 



2l{C>o}C 2 dC 
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cxp < — V2k1 



{3je{l,...,fe}, V27Ce(^)>Zj} 

max{0, s(xj)} 
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3=1 k 



This shows that the processes (HJ and (J2J have the same finite dimensional distributions. 

Let Co be the space of real valued continuous functions on X C M. d and $ = {fi}i>\ a Poisson point 
process on Co such that 

ifii(x) = V2ir(iSi(x), i = l,2,..., 

with {Ci}i>i and £j are defined as in (fTJ). To shorten the notations we write /(x) = {f(x\), . . . , f(xk)} 
for all (random) function /: X — > R and x = (xi, . . . , Xk) G X k . 

For x G X k , the Poisson point process {(fii(x-)}i>i defined on M. k has intensity measure 



/•OO 

A x (A) = / Pr[\/27rCe(x) G A]C _2 dC, A c M. k Borel set, 
Jo 

where e(x) is a centered multivariate normal random vector with covariance matrix E x = {p(x 
i, j = 1, . . . , k. Provided the covariance matrix E x is positive definite, it can be shown — sec Appendix fAl 
that the intensity measure A x is absolutely continuous with respect to the Lcbcsgue measure and that 
the density of A x is 



/)} 



A x (z) = 7r-( fe - 1 )/ 2 |E x r 1/2 ax(z)- (fc+1)/2 r 



fc + 1 
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where a x (z) = z T S x 1 z. 

The explic i t expr ession of the density A x (z) is the starting point of the methodology developed in 
Dombrv et all |201lj where it is shown that the quotient 



A s | x ,z( u ) 



A( s . x )(u, z) 
A x (z) ' 



u e 



characterizes the density of the extremal functions, i.e., the atoms of $ that reach (at least) one hitting 
bound Z(xi) = Zi for some i £ {1, . . . , fc}. Consequently to be able to draw conditional simulations from 
the process (TT]), the distributions of the extremal functions have to be identified. 

Proposition 1. For (s,x) £ X m+k , z £ IR fe and provided that the covariance matrix £( SjX ) of the 
random vector £:{(s,x)} is positive definite, the function A s | XjZ (u),, u £ R m , corresponds to the density 
of a multivariate Student distribution with k + 1 degrees of freedom, location parameter 



— S S:X S X 1 z, 



and scale matrix 



~ a x (z) ,„ 



k + 1 

Proof. For all u £ R m we have 

A s |x,z(u) - 7T | Sx |_ 1/2 
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We start by focusing on the ratio a( sx )(u,z)/a x (z). Since 









S x:s 







(S s — S S:X S X 1 E X:S ) 1 (S s S S:X I] X 1 E X:S ) ^^Sg^S,;. 1 
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straightforward algebra shows that 



a (s,x)(u,z) 
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(u-^S-^u-zi) 



A* — E s:x S x z, 



Ox(z) " ' fc+1 ' " _= ^ _X ~ fc+1 

We now try to simplify the ratio |£( SiX )|/|£ x |. Using the fact that 



S — r S S:X S X 1 S X:S ) 
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where Id is the identity matrix, is a matrix with zero entries, combined with some more algebra yields 



J (s,x)| 



5js:xE x S x:s | 



fc + 1 

a x (z) 



Using the two previous results it is easily found that 

(u-n) T Z-Hu-ri 



-(m+fe+l)/2 



-fc+1" 



A s , x , z (u) = 7T- m ^(k + l)-Wa|s|-Va J i + 



fc + 1 



r( Mi ; 



which corresponds to the density of a multivariate Student distribution with the expected parameters. □ 

At first sight it seems intriguing that the extremal functions involve a multivariate Student distribu- 
tion. However this result is not as surprising as it might be at first glance. Indeed consider the following 
simple max-stablc process 



1 /2 

Z{x) = cv 27rmax(iXj; max{0, £,(a;)}, 

Z>1 



E 



X 1 ' 2 



v > 1, 
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where the normalizing constant c appears to ensure unit Frechet margins, {Ci}i>i are the points of a 
Poisson process on (0, oo) with intensity dA(£) = £~ 2 d£, Xi are independent copies of a random variable 
X such that v/X follows a x 2 distribution with v degrees of freedom and e% are independent replicates of 
a standard Gaussian process as in ([1}. All these random objects are assumed to be mutually independent. 
It is easily shown that 



Pr[Z(xi) < zi, . . . , Z(x k ) < Zfc] = cxp 



-E 



- ^ xl/2 max{0 l£ (x,)} 



cxp 
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max 
j=i,...,k 



max{0, e(xj)} 



exp |-V2ttE 
Pr[Z(a;i) < z x 



max ■ 

3=1,-,* 



c{0,e(xj)} 



,Z{x k ) < z k ] 



and the original Schlather process defined through standard Gaussian processes bv ISchlather [2002 [ has 
an equivalent construction based on Student processes. 

A similar result as Proposition [T] holds for the Brown- Rcsnick process where A s | XjZ (u) corresponds to 
the density of a multivariate log-normal distribution instead of a multivariate Student distribution. In 
spite of this d i fferen ce, the method for conditional simulation of Brown-Rcsnick processes proposed by 
Dombrv et al. I |201lj can be adapted straightforwardly to obtain conditional simulations of a Schlather 
process. For consistency we briefly recall the procedure: 

Step 1 Draw a random partition 9 of the set {x\ 1 . . . , Xk} from the following distribution 

M 

Pr[0 = r | Z(x) = z] cx TT A x (z T ) / A x , x . K)dUj, 

J-i J i" '■ : 3 

where r is a partition of {x\, . . . , Xfc}, \t\ is the size of the partition, v T . = (i^: i G Ij) and 
v T ? = : i ^ Ij) with Ij = {i : x,i E Tj} for all fc-upple v living in X k or M. k ; 

Step 2 Given r = (n, . . . ,r^), draw £ independent extremal functions (^+, . 
sional distributions are 



whose finite dimen- 



Pr[(^+(s)edv|Z(x)=z,6» = T]cx |y l {u<ZT , } A (S:XT , ) | Xxj , ZTj (v,u)du|dv, s <G X m , m E N, 



where !{.} is the indicator function and set Z + = max(^|, . 



Step 3 Independently from the two previous steps, draw a Schlather process Z conditioned to the 
constraint Z~(x) < z by thinning the Poisson point process i.e., 

Z~ = max{<p e ip(x) < z} . 



Then, Z = max{Z + ,Z~} follows the conditional distribution of Z given Z(x.) = z. We remind that 
the partition size I denotes the number of extremal functions needed to satisfy the conditional event 
Z(x) = z. The extremal function ipf associated to the component t , of r satisfies the c onstraints 
Lpj(xi) = Zi if Xi E Tj and Lpj(xi) < Zi if Xi (fc Tj — see Theorem 1 in Dombrv et al. [2011 for more 
details. Consequently in Step 2 a rejection algorithm has to be used to ensure that the sampled extremal 
functions meet these conditions. 

For practical purposes, Step 1 is the more challenging since it amounts to sample from a discrete 
distribution whose state space corresponds to all possible partitions of the set {xi, . . . ,Xk}- Since the 
number of possible partitions of a set with k elements corresponds to Bell numbers, it results in a 
combinatorial explosion even for a moderate number of conditioning locatio ns. For instance whe n k = 10 
there arc around 116000 possibilities. To bypass this computational burden. iDombrv et al. I |201l| propose 
the use of a Gibbs sampler in Step 1 and in the remainder of this paper, we will use their Gibbs sampler 
whenever k is too large, e.g., k > 10. 
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Table 1: Spatial dependence structures of Schlather processes based on an isotropic powered exponential 
correlation function with scale parameter X and shape parameter k. The correlation function parameters 
are set to ensure that 0(100) = 1.5. 

Sample path properties 
9i : Wiggly 8 2 : Smooth 6*3 : Very smooth 
~A 208 144 128 

k 0.5 1 1.5 
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Figure 1: Three (unconditional) realizations of a Schlather process with standard Gumbel margins 
and extremal coefficient functions 9i,02 and #3 — from left to right. The squares correspond to the 15 
conditioning values that will be used in the simulation study. The right panel shows the associated 
extremal coefficient functions where the black, red and green lines correspond respectively to 9±, 62 and 

03- 



3 Simulation Study 



In this section we check whether our algorithm is able to produce realistic conditional simulations of 
Schlather processes. To this aim we consider an isotropic powered exponential correlation function 



p(h) = exp <^ - 



h>0, A > 0, < K < 2, 



and three different spatial dependence configurations as summarized in Table [T] Figure [2] shows one 
(unconditional) realization for each sample path configuration as well as the corresponding extremal 
coefficient functions. These realizations will serve as a basis for our conditioning events. 

To check if our sampling procedure is accurate and given a single conditional event Z(x) = z, wc 
generate 1000 conditional realizations of a Schlather process with standard Gumbel margins and extremal 
coefficient function 61, 62 and $3. Figure [5] plots the 0.025, 0.5 and 0.975 pointwise sample quantilcs from 
these 1000 conditional simulations. As expected it can be seen that the sample paths used to get the 
conditional events lay most of the time in the pointwise confidence intervals. In addition wc can see 
that the sample quantile paths inherit the regularity driven by the shape parameter k of the powered 
exponential correlation function and t hat there is les s varia bility in regions close to some conditioning 
locations. Contrary to the results of Dombrv et al.l (2011 for Brown-Resnick processes, the sample 
quantiles do not converge to that of a unit Gumbel distribution since the Schlather process is not ergodic 
and cannot reach independence as we go far away from any conditio ning location. 

Since the conditional process Z \ Z(x) = z is not max-stable Dombrv et all 20 11 . one need to 
integrate out with respect to the conditional event to retrieve the max-stability property. To this aim 
we generate 1000 independent conditional events Z(x) = z, x being fixed, and for each single condi- 
tional event one conditional realization of a S chlather process. F igure [3] compares the pairwisc extremal 



coefficient estimates using the F-madogram [Coolev et all 120061 ] to the theoretical extremal coefficient 
functions 9i,02, #3 with k = 5, 10, 15 conditioning locations. As expected, whatever the number of con- 
ditioning locations and the spatial dependence structures are, the (binned) pairwise estimates match the 
theoretical curves indicating that our conditional simulations have the right spatial dependence structure. 
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Figure 2: Pointwisc sample quantilcs estimated from 1000 conditional simulations of Schlather processes 
with standard Gumbcl margins and extremal coefficient function #i,#2 and 62 (left to right) and with 
k = 5, 10, 15 conditioning locations — top to bottom. The solid black lines show the pointwise 0.025, 0.5, 
0.975 sample quantilcs and the dashed grey lines that of a standard Gumbcl distribution. The squares 
show the conditional points {(#,, Zi)}i=i,...,fc and the solid grey lines correspond to the simulated paths 
of Figure QJ 
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Figure 3: Comparison of the pairwise extremal coefheient estimates (using a binned F-madogram with 
250 bins on 1000 independent conditional simulations) and the theoretical extremal coefficient functions 
with k = 5, 10 and 15 conditioning locations — from left to right. The 'o', and 'x' symbols correspond 
respectively to the spatial dependence conhgurations 8± , 82 and 0% . The solid grey lines correspond to 
the theoretical extremal coefficient functions. 



G 




Figure 4: Left: Topographical map of Switzerland showing the sites and altitudes in metres above sea 
level of 16 weather stations for which annual maxima temperature data are available. Middle: Times 
series of the daily maxima temperatures at the 16 weather stations for year 2003. The black squares, 
red circles and green triangles indicate the annual maxima that occurred respectively the 11th, 12th and 
13th of August. Right: Comparison between the fitted extremal coefficient function from a Schlather 
process (solid red line) and the pairwise extremal coefficient estimates (gray circles). The black circles 
denote binned estimates with 16 bins. 

Table 2: Empirical distribution of the partition size for the temperature data estimated from 10 parallel 
Markov chains of length 1 000. 



Partition size 1 2 3 4 5-16 

Empirical probabilities (%) 2.47 21.55 64.63 10.74 0.61 



4 Application 

In this section we apply our results to get co nditional simulations of extreme t emperatures. The data 
considered here were previously analyzed by iDavison and Gholamrezaed [2011 1 and consist in annual 



maximum temperatures r ecorded at 16 sites in Switze r land d uring the period 1961-2005 — see Figure @] 
Following the work of IDavison and Gholamrezaee [201lj . we consider a Schlather process with an 



isotropic powered exponential correlation function. Due to the lack of covariates others than longi- 
tude, latitude and elevation, we consider the following trend surfaces for the generalized extreme value 
distribution parameters 

fx(x) = + Pi,»a±t(x), (3) 

*{X) = /V: (4) 

= A),c + Pi,(alt(x), (5) 

where alt (a;) denotes the altitude above mean sea level (in kilometres) and {[i(x),a(x),£(x)} are the 
location, scale and shape parameters of the generalized extreme value distribution at location x. 

A preliminary study showed that the simultaneous fit of the marginal and dependence parameters 
induced a bias in the estimation of the spatial dependence structure probably owing to the use of too 
simple trend surfaces and we decided to fit separately the marginal parameters /3 V and the spatial 
dependence parameters A, k — the latter were estimated by maximizing the pairwise likelihood obtained 
by transforming the data to unit Frechet margins using the empirical distribution function. The spatial 
dependence parameter estimates are A = 260 (149) and k = 0.52 (0.12) and the corresponding fitted 
extremal coefficient function, which is similar to some extent to our test case #3 in Section [3J is shown 
in the right panel of Figure 

In year 2003, western Europe was hit by a severe heat wave believed to be the hottest one ever recorded 
since at most 1540 ("2003 European heat wave", Wikipedia: The Free Encyclopedia). Switzerland was 
also largely impacted by this severe extreme event since the nation wide record temperature of 41.5°C 
was recorded that year in Grono, Graubunden — near Lugano. Consequently for our analysis we use as 
the conditional event the maxima temperatures observed in summer 2003 — se e Figure |U Based on our 
fitted max-stable model and using the Gibbs sampler of Dombrv et al.l 2011 1, we simulate 10 parallel 
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Figure 5: From left to right: maps on a 64 x 64 grid of the pointwisc 0.025, 0.5 and 0.975 sample quantilcs 
for temperature (°C) obtained from 10000 conditional simulations of a Schlather process. The squares 
show the conditional locations and the conditional values. The right panel shows the temperatures 
anomalies, i.e., the difference between the pointwisc conditional medians and the pointwise unconditional 
medians estimated from the fitted trend surfaces ©-([5]). 



Markov chains of effective length 1000 — with a burn-in period of length 500 and a thinning lag of 100 
iterations. The empirical distribution of the partition size estimated from these Markov chains is shown 
in Table [2] We can see that around 90% of the time the conditional realizations were a consequence of at 
most three extremal functions. Since our original observations were not summer maxima but maximum 
daily values, a close inspection of the times series in year 2003 reveals that the hottest temperatures 
occurred between the 11th and 13th of August, see Figure 01 and, to some extent, corroborates the 
empirical distribution of Table [5] 

Figure [5] shows the 0.025,0.5 and 0.975 pointwise sample quantiles obtained from 10000 conditional 
simulations on a 64 x 64 grid. The conditional pointwise median provides an estimation of the temperature 
at a given location while the 0.025 and 0.975 pointwise sample quantiles provide 95% pointwise confidence 
intervals. As expected, we can see that the largest temperatures occurred in the plateau region of 
Switzerland while temperatures were appreciably cooler in the Alps. The right panel of Figure[5]shows the 
difference between the pointwise conditional medians and the unconditional pointwise medians estimated 
from our fitted trend surfaces ([S])-©. The differences range between 2.5°C and 4.75°C and, as expected, 
the largest differences occur in the plateau region of Switzerland. 



5 Conclusion 



In this paper we deri ved the (regular) con ditional distribution of the Schlather process and adapt the algo- 
rithm introduced bv lDombrv et aL [201 1 [ to get conditional realizations from this process. The proposed 
framework was tested on simulated data and an application on extreme temperatures in Switzerland was 
gi ven. Results show that our proc edur e gives accurat e simu lations. This work completes the one started 
bv lDombrv and Evi-Minkol 2011 1 and Dombrv et al. I [201l| and therefore enables the use of conditional 
simulations from wide ly used max-stab le models. Future works co uld focus on t he random set version of 
the Schlather process Schlather . 2002 as well as the Smith model [Smith . 1990j . Although the (regular) 
conditional distribution of the former could be found following the lines of this paper, trying to get closed 
forms for the latter is likely to be more challenging since the spectral measure of the Smith process is 
not regular. The algorithm s used for this work have been implemented in C within the R fr amework 
R Development Core TeamL 12011 and will be collected in the R package SpatialExtremes [Ribatetl . 
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A Disintegration of A x (v4) 

In this appendix we show that the intensity measure of the Poisson point process {y>i(x)}j>i defined on 
M. k is absolutely continuous with respect to the Lebesgue measure and give closed form for its density. 
It is straightforward to see that for all x <G X k and Borel set icM 1 

/>oo />oo /> 

Ax(A) - J Prh/2<E(x) G A]C 2 d( = J J l {v ^ CyeA} /x(y)dyC 2 dC, 

where / x denotes the density of the random vector e(x), i.e., a centered Gaussian random vector with 
covariance matrix E x . Now the change of variable z = \/2ir^y gives 

= (2»)-'|E Ii r" 2 / C °/ j »p(- i; ^z''E;'z)r < * +2 Mzd< 
= (2 Ir )-*|S x r^^ °exp (-^Vs;'*) C'dCd. 

= / A x (z)dz, 

where A x (z) = ^-( fe - 1 )/ 2 |S x |- 1 /2 ax ( z )-(fe+i)/2 r { (fc + 1)/2 } and a x (z) = z T ^ x 1 z. 
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